#_test.py用于对原数据进行裁剪，并获取经纬度，DEM，Land/SeaMask等信息
#内置模块
import sys
import os
import time
import datetime
import json as js
import tkinter as tk
import tkinter.messagebox as tkmegbox               #消息提示框,未使用
from multiprocessing import Lock,Manager            #多进程锁和共享数据
from concurrent.futures import ThreadPoolExecutor   #多线程____线程池
from concurrent.futures import ProcessPoolExecutor  #多进程____进程池
import multiprocessing

#第三方模块
import cv2
import numpy as np
import matlab.engine
from osgeo import gdal                              #用于读取hdf等栅格结构数据
from netCDF4 import Dataset
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from matplotlib.font_manager import FontProperties
font=FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=14)

'''
https://www.mathworks.com/content/dam/mathworks/mathworks-dot-com/support/sysreq/files/python-compatibility.pdf
'''

#打开包含高程信息的*.hdf文件
#创建默认路径为当前文件路径的文件对话框,返回选取文件的绝对路径
def openhdf():
    from tkinter import filedialog
    windows=tk.Tk()
    f=filedialog.askopenfilename(filetypes=[('Contains elevation data','*.hdf'),('All Files','*.*')],
                                                   initialdir=sys.path[1])
    windows.destroy()
    if len(f)==0:
        print('用户已取消操作')
        input('回车键退出...')
        sys.exit(1)
    else:
        return f
    

#从打开的*.hdf文件获取指定区域高程数据
'''
lon_E     指定区域最东经度对应的数据索引值
lon_step  指定区域经度跨度对应的索引值跨度(由东向西)
lat_S     指定区域最北维度对应的数据索引值
lat_step  指定区域维度跨度对应的索引值跨度(由北向南)

经纬度转化为数据索引参见函数degree2index
'''
def gethight(lon_E,lon_step,lat_S,lat_step):
    hdf=gdal.Open(openhdf())
    #hdf=gdal.Open('test.hdf')
    #获取文件框架信息
    subdatasets=hdf.GetSubDatasets()
    """
    #本部分用于输出*.hdf文件的基本信息
    print('Number of subdatasets: {}'.format(len(subdatasets)))
    for sd in subdatasets:
    print('Name: {0}\nDescription:{1}\n'.format(*sd))
    """

    '''
        0: Longitude     [6001x6001] (32-bit floating-point)  (80,-160)
        1: Latitude      [6001x6001] (32-bit floating-point)  (60,-60)
        2: SensorZenith  [6001x6001] (32-bit integer)
        3: SensorAzimuth [6001x6001] (32-bit integer)
        4: Height        [6001x6001] (32-bit integer)
        5: Land/SeaMask  [6001x6001] (8-bit unsigned integer)
        6: LandCoverType [6001x6001] (8-bit unsigned integer)
    '''
    #读取指定框架(第4维)数据并切片
    #Wrong:*.hdf数据的坐标为(lat,lon),与*.nc数据坐标(lat,lon)相同;
    #但在裁剪时*.nc裁剪的函数使用的是(lon,lat);
    #为使得两个数据的经纬度一致,需要使用下面(lat,lon)的格式正确裁剪*.hdf数据
    #当然在数据进行索引时,均可使用(lat,lon)进行数据索引,详见Group1函数
    Height=(gdal.Open(subdatasets[4][0]).ReadAsArray())[lat_S:lat_S+lat_step,lon_E:lon_E+lon_step]
    Land_SeaMask=(gdal.Open(subdatasets[5][0]).ReadAsArray())[lat_S:lat_S+lat_step,lon_E:lon_E+lon_step]
    return Height,Land_SeaMask


#将输入的裁剪经纬度信息转化为数据索引信息
'''
#输入参数:
    longitude 经度信息
    latitude  维度信息
    lon_E     裁剪区域最东经度
    lon_step  裁剪区域经度跨度
    lat_S     裁剪区域最北维度
    lat_step  裁剪区域维度跨度
    eg. 裁剪区域110E-170W,30N-20S应当输入为110,80,30,50
#输出参数:
    lon       指定区域最东经度对应的数据索引值
    lon_step  指定区域经度跨度对应的索引值跨度
    lat       指定区域最北维度对应的数据索引值
    lat_step  指定区域维度跨度对应的索引值跨度

round()函数为四舍五入取整,当输入参数为小数时计算结果难免会产生偏差
'''
def degree2index(args):
    longitude=args[0][0]
    latitude=args[0][1]
    lon_E=args[1][0]
    lon_step=args[1][1]
    lat_S=args[1][2]
    lat_step=args[1][3]

    flag1=0;flag2=0
    longitude=np.array(longitude)
    latitude=np.array(latitude)
    for i in range(longitude.size):
        if flag1==0:
            if longitude[i]>=lon_E:
                lon=round(i)
                flag1=1
        elif flag2==0:
            if longitude[i]>=lon_E+lon_step:
                lon_step=round(i-lon+1)
                flag2=1
                
    for i in range(latitude.size):
        if flag1==1:
            if latitude[i]<=lat_S:
                lat=round(i)
                flag1=2
        elif flag2==1:
            if latitude[i]<=lat_S-lat_step:
                lat_step=round(i-lat+1)
                flag2=2
    #Wrong:正确的经纬度信息不需要要进行+1操作,由于所有研究数据均已裁剪完成,建议保持原状
    #args[1][0]=lon    正确经度
    args[1][0]=lon
    args[1][1]=lon_step
    #args[1][2]=lat    正确维度
    args[1][2]=lat
    args[1][3]=lat_step
    return args


#根据时间等参数创建文件列表
'''
#输入参数:
    old_rootpath  输入文件的根目录
    new_rootpath  保存文件的根目录
    start_time    文件开始日期 eg.'2020/1/9'
    end_time      文件结束日期 eg.'2021/2/18'
    day/night     第八时区时间 eg.day=5,night=15表示只获取当天北京时间凌晨5点到下午15点的数据
    
#输出参数:
    old_path      输入文件列表
    new_path      输出文件列表
'''
def creat_path(start_time,end_time,day,night):

    #将日期字符串解析为可读年月日
    start_time=start_time.split('/')
    start_time=[int(start_time[0]),int(start_time[1]),int(start_time[2])]
    end_time=end_time.split('/')
    end_time=[int(end_time[0]),int(end_time[1]),int(end_time[2])]
    
    #设置时间间隔,转化UTC时间并获得文件列表
    stepdays=10
    daymin=datetime.datetime(start_time[0],start_time[1],start_time[2],day)-datetime.timedelta(hours=8)
    daymax=datetime.datetime(start_time[0],start_time[1],start_time[2],night)-datetime.timedelta(hours=8)
    end_time=datetime.datetime(end_time[0],end_time[1],end_time[2],night)-datetime.timedelta(hours=8)
    filename=[]
    while daymin<=end_time:
        while daymin<=daymax:
            filename.append(daymin.strftime('%Y%m\\%d\\NC_H08_%Y%m%d_%H%M_R21_FLDK.06001_06001.nc'))
            daymin=daymin+datetime.timedelta(minutes=60)

            #daymin=daymin+datetime.timedelta(minutes=60)
            if daymin==daymax:
                filename.append(daymin.strftime('%Y%m\\%d\\NC_H08_%Y%m%d_%H%M_R21_FLDK.06001_06001.nc'))
                daymin=daymin+datetime.timedelta(hours=day-night+24*stepdays)
                daymax=daymax+datetime.timedelta(days=1*stepdays)
                break
    return filename


#调用matlab2018引擎裁剪数据
def cutnc_matlab2018(args):
    print('正在处理{}'.format(args[3].split('\\')[-1]))
    old_path=args[2][0]+args[3]
    new_path=args[2][1]+args[3]
    lon_E=float(args[1][0])
    lon_step=float(args[1][1])
    lat_S=float(args[1][2])
    lat_step=float(args[1][3])
    ed=args[-2]
    ing=args[-1]
            
    eng1=matlab.engine.start_matlab()
    eng1.nc2cut(str(old_path),str(new_path),lon_E,lon_step,lat_S,lat_step,nargout=0)#先经度后维度
    
    ed.append(args[3])
    for i in ed:
        ing.remove(i)

    lock=Lock()
    with lock:
        with open('.log\\downloading','w',encoding='utf-8') as f:
            js.dump(ing,f)
            f.close()
            
        with open('.log\\.log','w',encoding='utf-8') as f:
            f.write('已完成:\n')
            for i in ed:
                f.write('{}\n'.format(i.split('\\')[-1]))
            f.write('\n未完成:\n')
            for i in ing:
                f.write('{}\n'.format(i.split('\\')[-1]))
            f.close()
    
    print('{}已完成'.format(args[3].split('\\')[-1]))
    print('--'*50)
    eng1.quit()


#构建日志文件以及参数
def setargs():
    if not os.path.exists('.log\\.args'):
        if not os.path.exists('.log'):
            os.makedirs('.log')
            print('配置文件夹不存在,已为您自动创建...')

        nc=Dataset('NC_H08_20151129_0300_R21_FLDK.06001_06001.nc')
        longitude=nc.variables['longitude'][:]
        latitude=nc.variables['latitude'][:]
        nc.close()
        lon=[];lat=[]
        for i in longitude:
            lon.append(float(i))
        for j in latitude:
            lat.append(float(j))
            

        lon_E=109;lon_W=122
        lat_N=20;lat_S=43

        day=12;night=12
        rootpath=['A:\\jma\\netcdf\\','']
        times=[['2019/8/30','2019/8/30']]

        '''
        day=14;night=18
        rootpath=['A:\\jma\\netcdf\\','']
        times=[['2020/12/27','2020/12/27']]
        '''


        filename=[]
        for i in times:
            filenames=creat_path(i[0],i[1],day,night)
            for j in filenames:
                filename.append(j)

        args=[[lon,lat],                                #文件经纬度坐标
              [lon_E,lon_W-lon_E,lat_S,lat_S-lat_N],    #研究区域经纬度范围
              rootpath,                                 #文件路径
              filename]
        
        with open('.log\\.args','w',encoding='utf-8') as f:
            js.dump(args,f)
            f.close()
        print('已将配置保存于".log\\.args"')
        
        with open('.log\\downloading','w',encoding='utf-8') as f:
            temp=[]
            for i in args[3]:
                temp.append(i)
            js.dump(temp,f)
            f.close()
            
        
    else:
        with open('.log\\.args',encoding='utf-8') as f:
            args=js.load(f)
            f.close()
        print('已从".log\\.args"加载配置数据...\n')
        
    return args


#多线程检校数据裁剪的完整性
def mulThread(thread_num):
    path=args[2][1]
    filenames=args[3]
    size=[args[1][3],args[1][1]]
    
    start_time1=time.time()
    notpass=[]
    
    with ThreadPoolExecutor(thread_num) as Thread_executor:
        lens=len(filenames)
        j=0
        for i in filenames:
            j=j+1
            filepath=(path+i).split('.nc')[0]+'_cut.nc'
            arg=[filepath,size,lens,j]
            notpass=notpass+Thread_executor.submit(checkfiles,tuple(arg)).result()
    if notpass==[]:
        print('\n所有文件均已通过校检...')
    else:
        print('\n共有{}个文件未通过校检:'.format(len(notpass)))
        for i in notpass:
            print(i)
            
    end_time1=time.time()
    print('本次校检共耗时{}s'.format(end_time1-start_time1))


#多进程裁剪数据,结合matlab函数
def mulProcess(process_num,args):
    start_time1=time.time()
    manager=Manager()
    ed=manager.list();ing=manager.list()
    with open('.log\\downloading',encoding='utf-8') as f:
        ing=js.load(f)
        f.close()
    if ing==[]:
        print('所有文件均已处理完毕...')
        #sys.exit(0)
        
    args.append(ed)
    args.append(ing)

    with ProcessPoolExecutor(process_num) as Process_executor:
        for i in ing:
            args[3]=i
            Process_executor.submit(cutnc_matlab2018,tuple(args))
            
    end_time1=time.time()
    print('多进程耗时{}s\n'.format(end_time1-start_time1))

   
def checkfiles(arg):
    nc=Dataset(arg[0])
    testdata=nc.variables['Hour'][:][-1][-1]
    if str(type(testdata))=="<class 'numpy.ma.core.MaskedConstant'>":
        return [arg[0]]
    else:
        print('\r共计{}个文件,正在校验第{}个文件...'.format(arg[2],arg[3]),end='')
        return []


if __name__=='__main__':
    args=setargs()
    args=degree2index(args)
    
    print('CPU核数为{}'.format(multiprocessing.cpu_count()))
    #经测试,使用多进程处理较快
    mulProcess(3,args)#启动3个进行裁剪数据
    #mulThread(1)#启动1个线程检校数据

    '''
    if not os.path.exists(new_rootpath):
        os.makedirs(new_rootpath)
    '''
    #longitude,latitude,Hight,file_path=getargs()
    #Hight,Land_SeaMask=gethight(args[1][0],args[1][1],args[1][2],args[1][3])
    





        
